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ABSTRACT 


Progress  in  developing  a  tool  to  compute  large 
amplitude  ship  motions  is  reported.  In  particular,  a 
method  to  calculate  transient  two-dimensional  potential 
flow  about  a  body  moving  in  a  free  surface  is  described. 
The  flow  problem  is  formulated  as  an  initial -boundary 
value  problem  in  which  the  velocity  potential  along  the 
free  surface  and  the  positions  of  the  moving  boundaries 
are  sought  as  solutions  of  a  coupled  system  of 
differential  equations.  An  implicit  finite-difference 
method  is  used  to  advance  the  solution  of  the  coupled 
system  .of  equations  in  time.  The  auxiliary  problem  of 
computing  the  velocity  potential  inside  the  fluid  region 
is  solved  by  a  method  which  is  based  on  boundary-fitted 
coordinates  and  is  directly  extensible  to  three- 
dimensional  flows.  Results  from  calculating  the  potential 
flow  about  a  body  in  forced  heave  motion  are  presented. 
The  hydrodynamic  force  on  the  body  has  been  obtained  and 
compared  with  the  hydrodynamic  force  predicted  from 
second-order  perturbation  theory. 
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INTRODUCTION 

In  recent  years  attention  of  naval  architects  and  ocean  engineers  has 
focused  on  how  vesse  . s  and  offshore  structures  react  to  large-amplitude  ocean 
waves.  Th  attention  has  been  motivated  by  the  capsizing  of  vessels  in  large 
breaking  waves  and  structural  failure  due  to  the  slamming  forces  associated 
with  such  waves.  It  is  therefore  of  interest  to  have  a  method  available  to 


determii  •  the  forces  on  a  floating  body  and  how  the  body  will  react  in  these 
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extreme  conditions.  Since  little  is  known  of  how  a  ship  reacts  even  to  non¬ 
breaking  waves,  it  would  be  extremely  valuable  to  have  a  method  available 
for  predicting  ship  motions  in  the  presence  of  large-amplitude  nonbreaking 
waves.  The  method  would  be  of  practical  value  in  a  systematic  study  of  how 
ship  design  changes  would  add  stability  in  rough  seas. 

Researchers  have  spent  much  effort  in  devising  methods  for  computing 
large-amplitude  ship  motions.  Their  work  is  usually  based  on  the  assumptions 
that  the  fluid  is  incompressible  and  the  fluid  motion  is  irrotational .  The 
assumptions  lead  to  the  existence  of  a  velocity  potential,  which  simplifies 
the  problem  formulation.  But,  since  the  equations  describing  the  free  surface 
are  nonlinear  and  cannot  be  linearized  for  large-amplitude  waves,  great 
difficulties  arise  in  the  computation  of  solutions  to  free-surface  potential 
flow  problems.  When  a  body  is  present  in  the  f^ee  surface,  additional 
difficulties  related  to  the  intersection  of  the  free  surface  and  the  body 
occur.  For  instance,  the  potential  flow  in  the  region  is  known  to  be 
singular.  Lin  et  al .  [1]*  have  recently  described  some  aspects  of  the 
singularity.  Dagan  and  Tulin  [2]  and  Fernandez  [3)  discuss  nonl inear ' ties  in 
fluid  flow  about  blunt  bodies.  Because  of  the  formidable  difficulties,  most 
of  the  work  on  nonlinear  free  surface  flows  has  been  for  two-dimensional 
(2-D)  flows. 

Many  authors  have  formulated  the  2-D  problem  as  an  ini tial -boundary  value 
problem  whose  solution  is  obtained  from  an  integral  equation  for  functions 
defined  along  the  boundaries  of  the  fluid.  Longuet-Higgins  and  Cokelet  [4] 
used  such  a  method  combined  with  a  time-stepping  procedure  to  calculate  free- 
surface  heights  with  no  body  present.  Faltinsen  |5]  and  Vinje  and  Brevig  [6-8] 

*A  complete  listing  of  references  is  given  on  pages  39  and  40. 
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have  used  the  integral-equation  approach  in  their  studies  of  fluid  motion 
in  the  presence  of  bodies.  They  make  the  restrictive  assumption  that  the  fluid 
domain  is  periodic  and  use  complex-variable  techniques  that  cannot  be  extended 
to  three-dimensional  problems.  Greenhow  et  al .  [9]  have  applied  the  method 
of  Vinje  and  Brevig  to  the  capsizing  of  a  body  in  the  free  surface.  Baker  et 
al .  [10]  have  developed  a  generalized  vortex  integral-equation  technique  that 
has  been  used  for  a  body  under  the  free  surface.  It  is  not  clear  whether  the 
generalized  vortex  method  is  suitable  for  numerical  computations  when  a  body 
intersects  the  free  surface  or  whether  it  will  be  computationally  efficient 
when  it  is  extended  to  three  dimensions.  Thus,  even  for  computing  nonlinear 
two-dimensional  free-surface  potential  flows,  full  generality  has  not  been 
attained . 

This  report  describes  progress  made  in  developing  a  tool  to  compute 
large-amplitude  ship  motions.  The  method  discussed  is  based  on  an  initial¬ 
boundary  value  formulation,  and  it  is  a  method  that  is  directly  extensible  to 
three  dimensions.  The  velocity  potential  along  the  free  surface  and  the  posi¬ 
tions  of  the  moving  boundaries  are  sought  as  solutions  of  a  coupled  system  of 
differential  equations.  An  implicit  finite-difference  method  is  used  to  march 
the  solution  of  the  coupled  system  of  equations  forward  in  time.  The  auxiliary 
problem  of  computing  the  velocity  potential  inside  the  fluid  region  is  solved 
by  a  finite-difference  method  based  on  boundary-fitted  coordinates.  Haussling 
[11]  has  presented  a  review  of  such  techniques  used  for  fluid  flow  problems. 

Results  from  calculating  the  potential  flow  about  a  cylinder  in  forced 
heave  motion  are  presented.  The  hydrodynamic  force  on  the  body  has  been  ob¬ 
tained  and  compared  with  the  hydrodynamic  force  predicted  from  second-order 
perturbation  theory. 
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MATHEMATICAL  FORMULATION 

The  physical  flow  problem  is  to  compute  transient  two-dimensional  flow 
about  a  body  moving  in  a  free  surface.  It  is  formulated  mathematical ly  as  a 
potential  flow  problem  in  which  the  velocity  potential  along  the  free  surface 
and  the  position  of  the  moving  free  surface  are  sought  as  the  solution  to  a 
nonlinear  initial-boundary  value  problem. 

In  particular,  the  physical  problem  is  to  determine  the  fluid  motion 
caused  by  the  prescribed  movement  of  a  body  partially  submerged  in  a  fluid  and 
the  resulting  hydrodynamic  force  on  the  body.  The  prescribed  motions  are 
forced  harmonic  heave  motions  never  so  large  that  the  body  becomes  completely 
submerged  in  or  rises  out  of  the  fluid.  The  body  considered  in  this  report  is 
a  closed  U-shaped  cylinder.  Gravity  is  the  only  body  force  acting  on  the 
fluid  which  is  inviscid,  incompressible,  and  initially  at  rest.  The  fluid 
motion  is  irrotational  and  thus  a  velocity  potential  '  is  assumed  to  exist. 
Surface  tension  is  neglected. 

All  variables  are  nondimensionalized.  Lengths  are  scaled  by  a  length  L 
characterizing  the  size  of  the  body.  Time  is  scaled  by  1/c  where  a  is  the 
frequency  of  the  body  motion  in  radians  per  second.  Velocities  are  scaled  by 
aL;  the  velocity  potential  $ ' ,  by  aL^ ;  pressure,  by  pa ;  and  force,  by 
po^L^  .  Here  p  is  the  fluid  density.  Thus,  for  example: 


x  '  =  Lx, 


y'  =  Ly, 


t '  =  t/o  , 


d) '  =  o L^d  ,  p'  =  pa2p2pi  p'  =  pa^L^F, 

where  (x,y)  are  variables  representing  the  coordinate  system,  t  is  time,  p  is 


pressure,  and  F  is  force.  The  primed  variables  represent  dimensional 
quantities;  the  nonprimed  variables,  nondimensional  quantities. 

The  fluid  region  is  described  in  terms  of  a  fixed  (x.y)-coordinate  system 
chosen  so  that  the  y-axis  points  vertically  upward  and  the  undisturbed  free 
surface  is  at  y  =  0  (Fig.  1).  A  fluid  region  of  infinite  depth  and  infinite 
lateral  extent  is  modeled  by  a  rectangular  tank  so  deep  that  the  effect  of  the 
bottom  boundary  is  insignificant  and  so  wide  that  no  waves  reflect  from  the 
side  boundaries  during  the  time  for  which  the  fluid  motion  is  modeled.  The 
rectangular  tank  is  bounded  by  the  lines  y  =  -h,  x  =  w,  and  x  =  -w.  The 
contour  of  the  body  moving  in  the  free  surface  is  given  as  a  function  of  time 
t  and  a  parameter,  a,  by  the  equations 

x  =  Bx(a,t)  =  A  (cos  (a)  -  0.1  cos  (3  a))  (la) 

y  =  By ( a , t )  =  A  (sin  (a)  +  0.1  sin  (3  a))  -  h0  cos  (t)  (lb) 

where  A  is  a  measure  of  the  size  of  the  body,  a  is  the  angle  measured 
counterclockwise  from  the  direction  of  the  positive  x-axis  (Fig.  2),  and  h0  is 
the  amplitude  of  the  heave  motion.  The  position  of  the  free  surface  is  given 
in  terms  of  a  parameter  e  and  the  time  t  by  x  =  xp(e,t)  and  y  =  yp(e,t).  The 
functions  xp(e,t)  and  yp(e,t)  are  to  be  calculated. 

Since  the  flows  considered  in  this  paper  are  symmetric  about  the  y-axis, 
only  the  half  of  the  fluid  region  where  x  >  0  is  considered  (Fig.  3).  The 
region  is  bounded  by  five  curves.  Across  the  boundaries  AE  (x  =  0), 
CD  (x  =  w),  and  ED  (y  =  -h)  there  is  no  flow.  The  curved  line  AB,  given  by 
Equations  (la,b),  is  the  contour  of  the  moving  body.  BC  is  the  tree  surface, 


whose  location  must  be  computed. 


The  velocity  potential  satisfies  the  Laplace  equation 

=  0  (2) 

in  the  fluid  region  and  is  subject  to  certain  boundary  conditions.  (Fluid 
velocity  in  the  x-direction  is  given  by  u  =  <|>x;  fluid  velocity  in  the  y- 
direction,  by  v  =  <j>y.)  At  a  solid  boundary,  the  normal  velocity  of  the 
fluid  must  equal  the  normal  velocity  of  the  solid  boundary  since  fluid  cannot 
penetrate  the  boundary  and  no  cavities  are  assumed  to  form  in  the  fluid.  In 
particular,  at  stationary  boundaries  the  normal  velocity  must  vanish.  At  the 
right  vertical  boundary  of  the  flow  domain,  about  16  half-beams  away  from  the 
body,  this  condition  is  given  by 

4>x  =  0atx  =  w  (3) 

Similarly,  at  the  bottom  boundary,  about  six  half-beams  below  the  free 

surface,  vanishing  normal  velocity  is  specified  by  the  equation 

4>y  =  0  at  y  =  -h  (4) 

At  the  boundary  AE  directly  beneath  the  body,  where  a  symmetry  condition  is 
specified  as  a  wall  condition,  the  velocity  potential  must  satisfy  the 
equation 

4>x  =  0  at  x  =  0  (5) 

Along  the  body  contour,  the  normal  velocity  is  known  from  the  prescribed 
motion  of  the  body.  In  fact,  the  normal  velocity  of  the  body  at 

(BxU.t),  By(a , t) )  is 

vn(a,t)  =  nx3Bx/8t  +  ny3By/3t 
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where  "ni(a,t)  =  (nx,ny)  is  the  unit  normal  directed  into  the  fluid  given  by 

(nx,ny)  =  (3By/3a,  -3Bx/3a)/[(3Bx/3a)2  +  (3By/3a) 2] 1/2 

The  required  boundary  condition  for  <|>  at  (Bx(a,t),  By(a,t))  on  the  body 
contour  is  thus 

3<)> 

—  =  4>xnx  +  4*yny  =  nx3Bx/3t  +  ny3By/3t  (6) 

3  n 

The  free-surface  coordinates  xp(e,t)  and  yp(e,t)  have  been 

parameterized  in  terms  of  e.  The  parameterization  e  is  chosen  such  that  for 
fixed  e  the  functions  xp(e,t)  and  yp(e,t)  describe  the  path  of  a  fluid 

particle.  In  other  words,  e  is  a  Lagrangian  variable.  The  velocity  potential 
on  the  free  surface  is  also  parameterized  in  terms  of  the  Lagrangian  variable 
e  by  =  4>p(e,t).  An  equation  to  be  satisfied  by  <j>p(e,t)  can  be  obtained 
from  Bernoulli's  equation,  which  can  be  expressed  as 

p  +  3*/3t  +  (4>2  +  <j>2)/2  +  (g/o2L)y  =  0  (7) 

x  y 

where  p  is  pressure,  g  is  the  acceleration  of  gravity,  a  is  the  frequency  of 
the  forced  harmonic  heaving,  and  L  is  the  characteristic  length.  Bernoulli's 

equation  is  valid  on  the  free  surface  and  throughout  the  fluid  region.  At  the 

free  surface,  the  pressure  is  assumed  to  be  zero  and  a  particular  case  of 
Bernoulli's  equation,  the  dynamic  free-surface  boundary  condition,  results: 

— F  (e , t )  =  —  (xF(e,t),yF(e,t),t)  =  (<j>2  +  4>2 )/2  -  (g/o2L)  yF(e,t)  (8) 

dt  Dt  r  r  x  y 

Here  D/Dt  =  4>x3/3x  +  4>y3/3y  +  3/3t  is  the  derivative  following  the  motion 
of  a  fluid  particle.  The  kinematic  free-surface  boundary  condition,  which 


states  that  no  fluid  particle  on  the  free  surface  can  leave  the  free  surface, 
is  expressed  by  the  two  equations 
dx-,  Dx 

— L  (e, t)  =  —  (xF(e,t),  yF(e,t),  t)  =  <f>  (xF(e,t),  yF(e,t),  t)  (9) 
dt  Dt  r  r  x  r  f 

dyF  Dy 

— -  (e,t)  =  —  (xp(e ,  t) ,  yF(e,t) ,  t)  =  $  (xp(e,t),  yF(e,t),  t)  (10) 
dt  Dt  r  r  y 

At  t  =  0,  the  velocity  potential  on  the  free  surface  is  given  by 

<f>F(e,t=0)  =  0  (11) 

and  the  free  surface  is  such  that 

yF(e, t=0)  =  0  (12) 

The  parameter  e  and  the  function  xF(e,t)  can  be  arranged  so  that 

xF(e, t=0)  =  e  (13) 

It  is  convenient  to  parameterize  the  fluid  at  the  body  contour  by  a 
Lagrangian  variable  e  so  that  x  =  xjj(e,t)  and  y  =  yg(e,t)  along  this 
boundary.  This  is  possible  since  fluid  particles  along  the  solid  boundary  can 
never  leave  that  boundary,  except  possibly  to  become  free-surface  particles  at 
the  intersection  of  the  body  and  the  free  surface.  Thus 

dxR  Dx 

— “-(e.t)  =  —  (xB(e,t),  yB(e,t),  t)  =  <t>x(*B(e,  t) ,  yB(e,t),  t)  (14) 


dyj 

A  + 


Dy 


(e,t)  =  —  (xB(e,t),  yB(e,t),  t)  =  <t>  (xB(e,t),  yB(e,t),  t)  (15) 


subject  to  the  initial  conditions 


xg(e,0)  =  BxCaCe) ,0) 
yB(e,0)  =  By(a(e) ,0) 


where  a(e)  is  a  prescribed  function. 

In  summary,  we  seek  the  solution  of  an  initial-boundary  value  problem  for 
x  =  xp(e,t),  y  ■  yp(e,t),  =  4>p(e,t)  on  the  free  surface  and  x  =  xg(e,t), 
y  =  yg(e,t)  on  the  body  contour,  in  which  e  is  a  Lagrangian  variable  that 
parameterizes  the  free  surface  and  the  body  contour.  These  five  functions 
obey  the  evolution  Equations  (8),  (9),  (10),  (14),  and  (15)  subject  to  the 
initial  conditions  (11),  (12),  (13),  (16),  and  (17).  The  velocity  potential 
in  these  equations  must  satisfy  Equation  (2)  subject  to  the  Neumann  boundary 
conditions  (3),  (4),  (5),  and  (6)  and  a  Dirichlet  boundary  condition  along  the 
free  surface  governed  by  Equation  (8). 

The  force  on  the  body  is  calculated  by  integrating  the  pressure  over  the 
wetted  surface  of  the  body.  Because  the  flow  problem  is  symmetric  about 
x  =  0,  the  x-component  of  the  force  on  the  body  vanishes: 


Fx  =  0 


Ine  y-component  of  the  force,  positive  upward,  is  given  by 


-2  J*  p(a,t)  ny(a,t)  (ds/da)  da 


where  (nx  ,ny)  is  the  unit  normal  at  the  body  contour  directed  into  the  fluid 
and  s,  increasing  in  the  counterclockwise  direction  along  the  body,  represents 
arclength.  Because  of  symmetry,  the  pressure  is  integrated  over  the  half  the 


vCaOcCfti 


wetted  length  of  the  body  given  by  x  =  Bj^a.t)  and  y  =  By(a,t)  where 
-  ir/2  <  a  <  aR(t).  (The  function  ap(t)  depends  on  the  position  of  the 
intersection  of  the  free  surface  and  the  body  contour.) 


NUMERICAL  SCHEME 


The  functions  xp(e,t),  yp(e,t),  $p(e,t)  along  the  free  surface  and 
xg(e,t),  yg(e,t)  along  the  body  contour  obey  five  coupled  first  order 
differential  equations  in  time  with  specified  initial  conditions.  In 
addition,  the  velocity  potential  must  satisfy  specified  boundary  conditions  at 
all  times.  To  solve  these  equations  a  finite-difference  method  is  used. 


BOUNDARY  FUNCTIONS 


Each  of  the  five  boundary  functions  is  discretized  with  respect  to  time 
and  space  using  a  fixed  time  step  At.  The  discretized  forms  of  the  functions 
are  denoted  by 


for  j  =  I, 


.  ,  N,  and 


x<n>  =  xp(e,,  nAt) 
Fj  F  J 


=  yv(ei»  nAt> 

Fi  t  J 


nAt) 

Fj  f  J 


x<n>  =  xR(ek,  nAt) 
Bk  B  k 


(20a) 


(20b) 


(20c) 


( 20d ) 
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'v>M  ■ 


y ^  =  yB(ek«  nAt> 


( 20e) 


for  k  =  1,  .  .  .  ,  M.  (The  subscript  j  will  be  used  for  a  free  surface 
variable  and  the  subscript  k  for  a  variable  along  the  body  contour.)  Thus,  the 
boundary  functions  are  discretized  into  boundary  grid  points  that  are  now 
assumed  to  move  as  if  they  were  associated  with  fluid  particles.  The 
evolution  Equations  (8),  (9),  and  (10)  for  $p(e,t),  xp(e,t),  yp(e,t)  and  the 
evolution  Equations  (14)  and  (15)  for  xg(e,t),  yg(e,t)  are  applicable  to  these 
particles  and  are  replaced  by  finite  difference  equations  based  on  the  Euler- 
modified  method.  The  finite-difference  equations  are  given  by 


(n+1) 

=  x(h) 

+  At(4>(n>  +  4>Cn+l>  ]/2 

(21a) 

Fj 

Fj 

xj  xj 

(n+1) 

=  y(n) 

+  At  (4>(n)  +  <t>(n+D]/2 

(21b) 

Fj  Fj  yj  yj 


4,(n+l) 

Fj 


x 


y 


-  (g/Lo2) 

fy  <  n  )  . 

y(n+!  )l\  /2 

(21c) 

lF) 

FJ  v 

(n+1) 

=  x<n>  + 

At 

+  (n) 

+  jjCn+l  ) 

!  ') 

(  2  Id  ) 

Bk 

Bk 

xk 

xk 

(n  +  1) 

=  y(n)  + 

At 

( n ) 

+  r»+ 1 ) 

f  2 

(  2  1  v  ) 

Bk 

Bk 

yk 
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J 

where 


^(m)  =  ^  ( x  ( m 1 ,  y  ( m )  ,  mAt) 

xl  x  Cl  Cl 


(j>(m)  =  ^  (x(m) ,  y(m) ,  mAt) 
yf.  y  Cl  (A 


1  1 
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for  i  -  j  or  k,  m  ■  n  or  n+1,  and  G  *  B  or  F.  The  initial  conditions  become 


m 

'  -  n 

.<-*■ 

h*.v 


x(0)  *  e 

(22a) 

Fj  j 

y(0)  =  o 

(22b) 

Fj 

a>(0)  =  o 

(22c) 

Fj 

(0) 

=  x  (a(e  ),  0) 

(22d) 

Bk 

B  k 

(0) 

=  y  (a(e  ),  0) 

(22e) 

Bk 

B  k 

In  Equations  (21a-e)  and  Equations  (22a-e),  the  subscript  j  runs  over  all 
possible  free  surface  grid  points,  and  the  subscript  k  runs  over  all  possible 
body  grid  points.  The  intersection  of  the  free  surface  and  the  body  is 
treated  as  a  body  point  obeying  Equations  (21d,e)  subject  to  an  initial 
condition  given  by  Equations  (22d,e). 

The  numerical  scheme  is  implicit.  An  initial  estimate  for  the  five 
functions  at  the  (n+l)-st  time  step  is  obtained  by  linearly  extrapolating  from 
two  previous  time  steps.  (For  the  first  time  step,  the  initial  estimate  is 
the  initial  condition.)  The  functions  are  corrected  iteratively  by  using 
Equations  (21a-e).  The  iterative  procedure  is  stopped  when  the  x-  and  y- 
values  have  satisfied  an  absolute  error  criterion  of  the  form 

|  f  (n+1 ,1)  _  f(n+l,Ul)|  <  £i  (23) 

and  the  <t»-values  along  the  free  surface  have  satisfied  the  relative  error 
criterion  given  by 

|l  -  <j>(n+l»^)/<j)(n+l  ,£+l)  j  <  e ^  whenever  |  <f>(n+l » O  |  > 
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The  superscripts  (n+l,£)  and  (n+l,£+l)  in  Equations  (23)  and  (24)  refer  to  the 
£-th  and  (£+l)-st  corrected  solutions  of  the  system  of  differential  equations 


at  time  step  n+1.  Generally,  the  stopping  criteria  for  the  iterative 
procedure  in  the  Euler-modified  method  at  each  time  step  are  that  the  maximum 
absolute  change  in  the  x-  and  y-values  on  the  free  surface  and  the  body  is 
Ej  =  0.001  and  the  maximum  relative  change  in  the  <|>-values  on  the  free 
surface  is  £2  =  0.001  for  <)>-values  greater  in  absolute  value  than 
e  3  =  0.000001. 

Each  of  the  five  discrete  evolution  equations  has  <|>x  or  <py  on  the  right 
side.  In  particular,  to  compute  the  right  sides  of  Equations  (21a-e),  <|>x  and 
<f>y  along  the  moving  boundaries  of  the  fluid  region  must  be  computed  from  the 
solution  of  the  Laplace  equation  for  4>.  The  solution  method  for  this 
equation  is  described  in  the  next  section. 

To  prevent  numerical  instabilities  on  the  computed  free  surface  from 
arising,  a  linear  filtering  scheme  due  to  Shapiro  [12]  is  used.  The  filtering 
scheme  has  been  used  successfully  by  Ohring  and  Telste  [13],  Haussling  and 
Coleman  [14],  and  several  other  researchers.  It  has  been  applied  at  fixed 
intervals  of  time,  and  has  been  especially  helpful  near  the  intersection  of 
the  free  surface  with  the  body  contour. 

Unless  some  method  is  used  to  maintain  a  reasonable  grid  spacing  along 
the  free  surface  and  the  body  contour,  grid  points  will  congregate  in  some 
areas  and  become  sparse  in  other  areas.  A  method  is  used  to  keep  a  uniform 
distribution  of  points  along  the  body  contour  and  a  prescribed  distribution  of 
grid  points  along  the  free  surface.  The  prescribed  free-surface  distribution 
is  such  that  the  free-surface  length  between  the  body  and  the  j-th  free- 
surface  grid  point  is  a  constant  fraction  of  the  total  free-surface  length 
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between  the  body  and  the  outer  boundary.  The  scheme  allows  one  to  follow  the 
movement  of  boundary  grid  points  within  a  time  step  as  if  they  were  fluid 
particles  and  to  shift  the  grid  points  to  other  fluid  particles  at  the  end  of 


each  time  step.  To  shift  the  grid  points,  cubic  spline  interpolation  is  used 


to  fit  the  arclength  as  a  function  of  grid  point  number  along  the  free  surface 
and  the  body  contour.  The  positions  of  the  grid  points  along  the  free  surface 
and  the  hull  are  shifted  and  the  values  of  all  pertinent  functions 
interpolated,  using  the  cubic  splines,  to  their  new  values.  The 
redistribution  of  grid  points,  of  course,  affects  the  initial  guess  for  the 
Euler-modified  method  at  time  step  n+1.  However,  if  the  shifting  is  done 
every  time  step  and  the  time  step  is  sufficiently  small,  the  redistribution 
scheme  has  been  found  to  proceed  smoothly. 

Because  of  numerical  errors,  grid  points  cannot  be  expected  to  remain 
exactly  on  the  body  as  the  solution  of  the  initial -boundary  value  problem  is 
advanced  in  time.  Numerical  errors  arise  from  the  redistribution  scheme  and 
from  replacing  the  differential  equations  by  finite  difference  equations.  To 
correct  for  such  errors,  grid  points  that  move  off  the  body  are  shifted 
back  to  the  body.  This  is  accomplished  by  computing  the  counterclockwise 
angle  about  the  center  of  the  body  from  the  direction  of  the  positive  x-axis 
(Fig.  2).  Every  body  grid  point  at  a  location  slightly  off  the  body  surface 
at  a  certain  value  of  that  angle  is  relocated  to  the  point  on  the  body  having 
the  same  value  of  that  angle. 

At  the  intersection  of  the  free  surface  and  the  body  contour  difficulties 
arise.  At  this  point  both  the  free-surface  boundary  conditions  and  the 
boundary  conditions  associated  with  the  solid  boundary  apply  at  a  particular 


a  position  on  the  free  surface  or  may  move  to  a  position  on  the  body  below  the 


surface.  In  other  words,  the  intersection  point  may  not  move  with  the  fluid. 
Since  the  free-surface  and  hull  boundary  conditions  all  involve  time 
derivatives  following  the  fluid  motion,  they  cannot,  in  general,  be  directly 
applied  over  a  time  interval  to  predict  quantities  at  the  intersection  point. 
Special  methods  for  handling  this  point  must  be  developed.  However,  for 
heaving  motions  of  an  almost  wall-sided  body,  the  intersection  will  be 
essentially  a  fluid  particle.  Therefore  for  the  current  study,  the  body 
Equations  (14)  and  (15)  are  applied  directly  at  this  point.  Inaccuracies  in 
this  approach  will  become  apparent  in  the  form  of  a  deviation  from  zero  of  the 
pressure  at  the  intersection  point.  In  fact,  such  pressure  deviations  might 
be  used  in  a  method  to  more  accurately  follow  the  intersection  point  as  would 
be  necessary  for  more  complicated  body  shapes.  One  such  scheme  has  been 
tested  but  has  proved  to  be  numerically  unstable. 

The  pressure  at  all  the  grid  points  along  the  body  including  the 
intersection  is  calculated  from  a  finite-difference  version  of  Bernoulli's 
equation: 


The  pressure  is  the  pressure  at  the  k-th  fluid  particle  on  the  body  at  time 
t  =  (n+1/2)  At.  The  force  on  the  body  at  this  time  is  calculated  by 
numerically  integrating  the  pressure  along  the  body  contour.  Trapezoidal 


quadrature  is  used 


LAPLACE  SOLVER 


To  solve  the  Laplace  equation  for  a  finite  difference  method  based 
on  boundary  fitted  coordinates,  due  to  Thompson  et  al.  [15],  is  chosen.  The 
f inite-dif ference  method  involves  mapping  the  time-dependent  fluid  region  onto 
a  fixed  computational  region.  The  coordinates  £  and  n  in  the  computational 


region  are  such  that 

they  obey  the  Laplace 

equation  with  x  and  y  as 

dependent 

variables: 

^xx  +  ?yy 

=  0 

(26a) 

^XX  +  hyy 

=  0 

(26b) 

The  boundary  conditions  for  £  and  n  along  a  given  boundary  are  Dirichlet 
if  a  particular  mesh  distribution  along  the  boundary  is  prescribed.  The 
boundary  condition  of  one  coordinate  is  Dirichlet  and  that  of  the  other 
coordinate  is  Neumann  if  mesh  orthogonality  near  a  particular  boundary  is 
desired. 

Since  all  computations  are  to  be  done  in  the  fixed  (£ ,n )-computational 


£ 


region,  it  is  convenient  to  interchange  the  independent  and  dependent 
variables.  When  this  is  done,  the  Equations  (26a, b)  for  £  and  n  become 


-  26x  +  yx 

=  0 

(27a) 

5n 

nn 

-  26y  +  yy 

=  0 

(27b) 

Sn 

nn 

where 


a  =  x2  +  y2 


6  =  x  x  +  y  y 
5  n  £  n 


Y  =  x2  +  y2 
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(27c) 

(27d) 

(27e) 
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To  obtain  a  mesh  that  wraps  around  the  body  and  conforms  to  the  other 
boundaries,  the  physical  fluid  region  is  divided  into  several  subregions  (Fig. 
4).  Each  subregion  is  mapped  onto  a  rectangle  of  computational  space.  In 
each  rectangle  of  computational  space,  the  inverted  Laplace  equation  is  solved 
subject  to  Dirichlet  boundary  conditions,  transformed  Neumann  boundary 
conditions,  or  matching  boundary  conditions  where  rectangles  overlap. 

The  Laplace  equation  for  the  velocity  potential  <f>  transforms  exactly  as 
Equations  (26a, b).  The  transformed  Laplace  equation  is  given  by 

cn(>  -  2fSj>  +  y  4.  =0  ( 27 f ) 

SS  nn 

where  a,  6,  Y«  are  defined  by  Equations  (27c-e). 

Wherever  possible,  central  differencing  is  used  to  discretize  Equations 
(27a-f).  At  the  boundaries  of  the  fluid  region  where  a  Neumann  boundary 
condition  is  specified,  second-order  one-sided  finite  differences  replace  some 
of  the  derivatives  in  Equations  (27a-f).  (See  Coleman  and  Haussling  [16]  for 
more  details.)  The  resulting  system  of  quasilinear  equations  for  x,  y,  and  $ 
at  the  grid  points  is  solved  by  successive  overrelaxation. 

Lin  et  al .  [1]  cite  the  works  of  various  researchers  who  show  that  a 
logarithmic  singularity  exists  in  the  velocity  potential  of  free-surface 
potential  flow  near  the  intersection  of  the  free  surface  witli  a  vertical 
wavemaker  in  horizontal  motion.  But,  since  the  body  contour  for  the  problem 
considered  is  nearly  vertical  at  the  intersection  with  the  tree  surface  and 
since  the  horizontal  velocity  component  of  the  body  is  zero,  the  logarithmic 
singularity  in  the  velocity  potential  may  be  relatively  weak.  Thus  special 
numerical  treatment  of  the  singularity  may  not  be  critical.  in  tact,  nothing 
special  has  been  included  in  the  numerical  method  to  accommodate  such  a 


singularity  if  it  should  arise. 


RESULTS 


Several  forced  harmonic  heave  motions  of  the  U-shaped  body  in  the  free 
surface  have  been  considered.  The  shape  of  the  body  has  been  fixed  by  setting 
the  parameter  A  in  Equations  (la,b)  to  0.7407  in  all  cases.  With  A  set  to 
this  value,  the  half-beam  of  the  body  and  the  draft  of  the  body  are  both 
0.6667.  Amplitudes  considered  were  hg/b  =  0.05,  0.3,  and  0.4,  in  which  hg  is 
the  amplitude  of  the  motion  and  b  is  the  half-beam  of  the  body.  Most 
researchers  consider  values  of  the  frequency  parameter  bo2/g  that  lie  between 
0.0  and  2.0.  In  this  study  the  frequency  parameter  is  restricted  to  lie  in 
the  interval  from  1.5  to  2.0  since  this  interval  contains  the  frequencies  for 
which  nonlinear  effects  are  greater. 

Linear  theory  for  the  problem  of  an  oscillating  body  in  a  fluid  of 
infinite  depth  and  lateral  extent  predicts  that  the  wavelength  far  from  the 
body  will  approach 

A/L  =  2ng/Lo2  =  2u(b/L)g/ba2 

asymptotically  in  time  [17].  The  dimensions  of  the  rectangular  tank  are  taken 
to  be  about  one  such  wavelength  deep  and  four  wavelengths  in  half-length. 
Thus,  if  the  frequency  parameter  bo2/g  is  no  smaller  than  1,  the  depth  h 
should  be  about  4  and  the  length  should  be  about  16.  This  region  is  long 
enough  that  no  waves  will  reach  the  side  boundary  during  the  first  few  periods 
of  forced  motion.  Since  the  region  is  more  than  half  a  wavelength  deep,  the 
effects  of  finite  depth  can  be  ignored. 

The  fluid  region  has  been  divided  into  five  time-dependent  subregions. 
Each  subregion  has  been  mapped  onto  a  fixed  rectangle  of  computational  space 
(Figs.  4,  5).  The  size  of  the  mesh  in  each  of  these  subregions  has  been 
depicted  in  Figure  5.  A  total  of  2916  grid  points,  counting  twice  those 
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points  where  two  subregions  overlap,  has  been  used  for  the  grids  that  cover 
the  computational  and  fluid  regions.  There  are  34  equally  spaced  grid  points 
on  the  half-body  contour.  At  the  far  right  vertical  boundary,  at  the  bottom 
horizontal  boundary,  and  at  the  vertical  boundary  below  the  body,  grid  points 
position  themselves  in  a  way  that  makes  the  mesh  near  these  boundaries 
orthogonal.  At  the  far  right  boundary  there  are  25  grid  points;  at  the  bottom 
boundary,  102  grid  points;  on  the  vertical  boundary  directly  beneath  the 
body,  where  a  symmetry  condition  is  specified  as  a  wall  condition,  18  points. 
The  initial  distribution  of  grid  points  on  the  free  surface  is  arranged  so 
that  the  mesh  near  the  intersection  of  the  free  surface  and  the  body  is 
approximately  uniform  in  both  the  x-  and  y-directions  (Fig.  6),  and  the  mesh 
near  the  intersection  of  the  free  surface  and  the  far  right  boundary  is  also 
approximately  uniform  in  both  directions.  There  are  90  grid  points  along  the 
free  surface. 

The  iteration  for  the  solution  to  the  Laplace  equations  for  the  changing 
mesh  and  for  the  velocity  potential  is  done  simultaneously.  It  is  found  that 
the  mesh  does  not  change  greatly  from  time  step  to  time  step  and  hence  a  test 
for  the  convergence  of  <£  suffices.  This  convergence  test  is  that  the  square 
root  of  the  sum  of  the  squares  of  the  residuals  at  the  34  body  points  should 
be  less  than  0.001.  A  spot  check  of  the  solution  iterates  indicated  that  the 
relative  change  in  <j>  near  the  body  was  about  0.1  percent  when  the  criterion 
was  satisfied. 

For  the  solution  of  the  Laplace  equations,  an  overrelaxation  factor  of 
1.4  was  chosen  for  grid  points  inside  the  computational  region.  For  grid 
points  on  the  body,  where  a  Neumann  boundary  condition  is  specified,  such  a 
large  relaxation  factor  caused  instabilities  in  the  solution  process  for  <p . 


during  the  first  two  periods  of  forced  motion  with  the  amplitude  such  that 


hg/b  =  0.4  and  the  frequency  such  that  ba^/g  =  2.  A  time  step  of  At  =  0.02 
resolved  one  period  of  motion  into  about  314  time  steps.  Figures  7  and  8  show 
the  free  surface  and  the  body  with  a  fixed  coordinate  system.  Figure  7, 
corresponding  to  times  between  5.4  and  8.8,  shows  a  rising  body  and  a  sinking 
free  surface  near  the  body;  Figure  8,  for  times  between  9.0  and  11.4,  shows  a 
descending  body  and  a  rising  free  surface.  Figures  9  and  10  show  similar 
results,  but  the  coordinate  system  in  these  figures  is  fixed  to  the  heaving 
body.  The  time-dependent  details  of  the  free  surface  near  the  body  are 
clearer  in  this  frame  of  reference.  Figure  9  shows  a  rising  body  between  the 
times  6.2  and  9.0;  Figure  10,  a  descending  body  between  the  times  9.2  and 
12.2.  It  is  interesting  to  note  that  the  slope  of  the  free  surface  near  the 
body  is  largest  at  about  the  time  the  body  has  attained  its  maximum  height. 
Results  for  computing  the  vertical  force  on  the  heaving  U-shaped  body  are 
depicted  in  Figure  11.  From  the  figure  it  is  seen  that  the  force  has  become 
periodic  in  less  than  one  period  of  forced  motion.  Also  shown  on  the  figure 
is  a  curve  of  the  force  versus  time  predicted  from  the  second-order  theory  of 
Papanikolaou  and  Nowacki  [18].  The  other  curve  is  a  linear  magnification  of 
the  results  from  calculating  the  force  when  the  amplitude  of  motion  is  eight 


times  smaller  (ho/b  =  0.05)  but  the  frequency  of  the  motion  is  the  same. 


Qualitative  agreement  among  the  curves  appears  to  be  good.  The  deviation 
between  the  computed  nonlinear  force  and  either  the  linear  or  the  second-order 
force  represents  aspects  of  force  which  could  not  be  predicted  from  linear  or 
second-order  theory. 

Another  comparison  of  the  computed  force  with  the  force  predicted  by 
second-order  theory  is  obtained  from  a  Fourier  analysis  of  the  curve  of  force 
versus  time  for  the  last  period  of  forced  motion.  According  to  second-order 
theory,  the  total  dimensional  force  on  the  body  when  its  center  (xg,yg) 
oscillates  vertically  according  to  the  formula  y  =  hg  sin  (at)  is  given  by 

F  =  2  pg  A  (0.97  tt/4)  +  2  pg  b^  e  Fj  sin  (at  +  6j) 

+  2  pg  b2  e2  F2Q  +  2  pg  b2  e2  F21  sin  (2at  +  62)  (28) 

where  e  =  hg/b  is  the  perturbation  parameter.  The  first  term  represents  the 
hydrostatic  force  on  the  body  at  its  neutral  water  position.  The  second  term 
represents  the  first-order  force  on  the  body,  and  the  last  two  terms  represent 
second-order  modifications  to  the  first-order  force  on  the  body.  A  Fourier 


analysis  of 

the  computed 

nonlinear 

force  versus 

time 

curve  will  produce 

coef  f i ci ent s 

of  the 

Fourier 

expansion 

in 

the 

orthogonal 

t unct ions 

{ 1 ,  sin  (nt ) 

,  cos  (nt)  }. 

The  coefficients 

are 

obtained  from 

the  usual 

integrals  by 

numerical 

integration 

.  In  the 

case 

of  the 

computed 

nonlinear 

force,  the  center  of  the  body  (xg,yg)  oscillates  vertically  according  to  the 
formula  yg  =  hg  cos  (ot)  instead  of  y  =  hg  sin  (at).  After  considering  this 
difference,  the  nonlinear  results  can  be  used  to  arrive  at  the  five  parameters 
describing  the  force  in  Equation  (28).  The  procedure,  of  course,  assumes  that 
the  contribution  to  the  force  from  third-order  effects  is  insignificant. 
Table  1  shows  how  the  computed  results  compare  with  those  of  Lee  [17],  and 


Papanikolaou  and  Nowacki  [18].  The  computed  results  seem  to  be  good  except 
for  the  second-order  amplitude  F2i«  The  discrepancy  In  the  second-order  force 
amplitude  is  probably  due  to  insufficient  accuracy  in  the  computed  nonlinear 
force.  Since  first-order  contributions  to  the  force  are  dominant,  the 

TABLE  1  -  FIRST-  AND  SECOND-ORDER  FORCE  COEFFICIENTS 


Papanikolaou 

and 

Nowacki 

Computed 

Nonlinear 

Results 

0.68 

0.67 

357° 

357° 

-0.19 

-0.23 

0.68 

0.82 

O 

O 

•—A 

108° 

relative  accuracy  of  the  remaining  force  when  the  first-order  predictions 
are  subtracted  is  small.  The  error  is  compounded  since  values  calculated  from 
the  remainder  are  divided  by  the  square  of  the  perturbation  parameter,  a 
small  number,  to  get  F20  and  F2 1 .  The  computed  second-order  phase  angle 


is.  however,  excellent. 


A  rough  estimate  of  the  accuracy  of  the  force  computations  can  be 
obtained  from  two  sets  of  computations  in  which  only  the  time  step  is 
different.  One  such  numerical  test  was  conducted  for  this  problem  with  the 
frequency  parameter  o^b/g  =1.7  and  the  amplitude  of  the  motion  such  that 
hg/b  =  0.3.  The  two  time  steps  used  were  At  =  0.02  and  At  =  0.01.  After 
about  t  =  2,  the  average  difference  in  the  computed  force  was  less  than  one 
percent.  Thus  the  error  seems  to  be  quite  small. 

Various  frequencies  of  forced  motion  have  been  considered  for  the 
amplitude  corresponding  to  hg/b  =  0.3.  Predictions  of  the  coefficients  for 
second-order  theory  were  made  and  compared  with  those  of  Papanikolaou  and 
Nowacki  [18],  Potash  [19],  and  Lee  [17].  The  results  are  presented  in  Figures 
12  and  13.  The  computed  first-order  coefficients  agree  well  with  the 
previously  computed  results.  The  computed  second-order  phase  angle  agrees 
well  with  the  results  of  Papanikolaou  and  Nowacki.  The  magnitude  of  the 
sinkage  (the  third  term  in  Equation  (28))  agrees  well  with  the  previously 
computed  results,  but  the  magnitude  of  the  oscillatory  part  of  the  second- 
order  force  does  not  agree  so  well.  This  is  a  reflection  of  the  relative 
error  in  the  nonlinear  force  when  the  first-order  force,  accounting  for  most 
of  the  force,  has  been  subtracted.  It  is  doubtful  that  there  is  more  than 
about  one  digit  of  accuracy  in  the  computed  second-order  results.  To  obtain 
these  results,  the  fluid  motion  resulting  from  about  two  periods  of  forced 
harmonic  oscillation  was  computed,  and  the  force  for  the  last  period  of  the 
motion  was  decomposed  into  its  Fourier  components. 
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CONCLUSION 


A  computer  program  to  compute  the  fluid  motion  resulting  from  forced 
harmonic  heaving  of  a  U-shaped  body  in  a  free  surface  has  been  produced.  The 
program  computes  forces  in  reasonable  agreement  with  the  results  of  previous 
researchers  who  used  second-order  perturbation  theory.  The  forces  have  been 
computed  for  cases  of  significant  nonlinear  fluid  motion. 

The  computer  program  will  be  modified  to  handle  a  variety  of  body  shapes 
and  motions.  Unlike  the  complex-variable  techniques  of  Faltinsen  [5]  and 
Vinje  and  Brevig  [6-8],  the  method  chosen  to  solve  the  problem  is  extensible 
to  the  problem  of  3-D  nonlinear  ship  motions.  The  fluid  domain  has  not  been 
assumed  periodic,  as  was  assumed  by  Vinje  and  Brevig  [6-8].  The  location  of 
the  intersection  of  the  free  surface  and  the  body  has  not  been  found  from 
extrapolation,  as  was  used  by  Vinje  and  Brevig  [7]. 

Additional  research  needs  to  be  performed  before  this  method  can  be 
applied  to  the  general  problem  of  ship  motions  in  two  or  three  dimensions. 
Satisfactory  results  have  been  obtained  for  a  heaving  U-shaped  cylinder,  but 
heaving  wedge-shaped  bodies  or  blunt  bodies  in  horizontal  motion  in  the  free 
surface  are  known  to  have  more  singular  fluid  flow  behavior.  Before  such 
behavior  can  be  treated  properly,  a  more  sophisticated  treatment  of  the  flow 
near  the  hull/water  surface  intersection  will  be  required.  Such  a  treatment 
should  include  both  the  dynamic  free-surface  equation  and  the  kinematic 
condition  associated  with  the  solid  boundary  of  the  cylinder. 
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X0  =  ° 

Vq  =  _hO  cos  W 


{x0  +  0.9A,  yQ) 


x  =  x0  +  A(cos  (a)  -  0.1  cos  (3a)) 
y  =  Vq  +  A(sin  (a)  +  0.1  sin  (3a)) 


Figure  2  -  Geometric  Description  of  the  Moving  Body  Contour 
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Surface  Positions  in  a  Frame  of  Reference  Fixed  to  the  Body 
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